setwd("E:/5hmc_file/motifbreakR_predict_enrich")

case=read.table("E:/5hmc_file/2_5hmc_yjp_bam/ASM/bayes_p/bias_AShM_BF_no_motif.txt",head=T,sep="\t")
case=case[case$BF_in_DC>10,]
#case=tidyr::separate_rows(case,Gene.refGene,sep=";")




con.file=list.files(path = "E:/5hmc_file/组织特异性表达/",pattern=".csv")
con.file=con.file[-c(grep("tissue",con.file))]
con.file=paste0("E:/5hmc_file/组织特异性表达/",con.file)

asm=data.frame(chr=case$chr,start=case$pos,ref=case$Ref,alt=case$Alt,rsid=case$avsnp150)
#snps=data.frame(chr=files2$Chr,start=files2$Start,ref=files2$Ref,alt=files2$Alt,rsid=files2$avsnp138)
asm$chr=gsub("chr","",asm$chr)
#snps$chr=gsub("chr","",snps$chr)




filekp=read.table("bias_AShM_BF_motif.txt",head=T,sep="\t")
 filekp=filekp[!filekp$avsnp150==".",]
filekp$seqnames=gsub("chr","",filekp$seqnames)
filekp$motif_region=paste(filekp$seqnames,filekp$start,filekp$end,sep=":")
symbol=unique(filekp$geneSymbol)
symbol=symbol[-c(1)]								#去掉第一个的NA


in_region=function(data_frame,chr,start,end){
part=data_frame[data_frame$chr==chr,]
parts=part[part$start<=end & part$start>=start,]
return(dim(parts)[1])
}
flanking_region=function(data_frame,chr,start,end){
part=data_frame[data_frame$chr==chr,]
parts=part[(part$start>end & part$start<=(end+500)) | (part$start<start & part$start>=(start-500)),]
return(dim(parts)[1])
}

result=data.frame(matrix(NA,1,ncol=8))
col_name=c("TFname","con.file","motif_ASM","motif_non_ASM","no_motif_ASM","no_motif_non_ASM","OR","p.value")
names(result)=col_name
for(k in con.file){
con=read.table(k,head=T,sep=",")
snps=data.frame(chr=con$Chr,start=con$Start,ref=con$Ref,alt=con$Alt,rsid=con$avsnp150)
snps$chr=gsub("chr","",snps$chr)
for(i in 1:length(symbol)){
	part =filekp[filekp$geneSymbol==symbol[i],]
	part=part[!is.na(part$unitID),]
	part_regions =unique(part$motif_region)
	motif_non_ASM=0
	motif_ASM=0
	no_motif_ASM=0
	no_motif_non_ASM=0
	for(j in part_regions){
		region=unlist(strsplit(j,":"))
		chr=as.numeric(region[1])
		start=as.numeric(region[2])
		end=as.numeric(region[3])
		a=in_region(asm,chr,start,end)
		b=in_region(snps,chr,start,end)
		c=flanking_region(asm,chr,start,end)
		d=flanking_region(snps,chr,start,end)
		motif_ASM=a+motif_ASM
		motif_non_ASM=b+motif_non_ASM
		no_motif_ASM=c+no_motif_ASM
		no_motif_non_ASM=d+no_motif_non_ASM
		}
		rt_tmp=data.frame(matrix(NA,1,ncol=8))
		names(rt_tmp)=col_name
		rt_tmp[,1] =symbol[i]
		rt_tmp[,2] =gsub(".hg19_multianno.csv","",gsub("E:/5hmc_file/组织特异性表达/","",k))
	rt_tmp[,3]=motif_ASM
	rt_tmp[,4]=motif_non_ASM
	rt_tmp[,5]=no_motif_ASM
	rt_tmp[,6]=no_motif_non_ASM
	rt_tmp[,7]=fisher.test(matrix(c(rt_tmp[,3],rt_tmp[,4],rt_tmp[,5],rt_tmp[,6]),nrow=2))$estimate
	rt_tmp[,8]=fisher.test(matrix(c(rt_tmp[,3],rt_tmp[,4],rt_tmp[,5],rt_tmp[,6]),nrow=2))$p.value
	result=rbind(result,rt_tmp)
	}
	}
	result$FDR =p.adjust(result$p.value,method = "BH")
	write.csv(result[-1,],"enrichment_by_TF_predit_by_MBR.csv",quote=F,row.names = F)
	
	
	
	######################################接下来是画图
	library(ggplot2)

file=read.csv("enrichment_by_TF_predit_by_MBR.csv",head=T)
#file[file$no_motif_AShM==0,]$OR=26
names(file)=c(names(file)[1:2],"counts",names(file)[4:9])
#file=file[order(file$OR),]
file1=file[file$p.value<0.05&file$OR>1,]
file2=file[file$FDR<0.1&file$OR>1,]
cf=unique(file$con.file)

file1$`1/OR`=1/file1$OR
file11=file1[file1$con.file==cf[1],]
file12=file1[file1$con.file==cf[2],]
file13=file1[file1$con.file==cf[3],]
file14=file1[file1$con.file==cf[4],]

#file11=file11[order(file11$p.value,decreasing = T),]
file11=file11[order(file11$OR),]
file11$num=1:dim(file11)[1]
p1=ggplot(file11,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+
  scale_y_continuous(breaks = file11$num,labels = file11$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())

#file12=file12[order(file12$p.value,decreasing = T),]
file12=file12[order(file12$OR),]
file12$num=1:dim(file12)[1]
p2=ggplot(file12,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+
  scale_y_continuous(breaks = file12$num,labels = file12$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p2

file13=file13[order(file13$OR),]
file13$num=1:dim(file13)[1]
p3=ggplot(file13,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+
  scale_y_continuous(breaks = file13$num,labels = file13$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p3

file14=file14[order(file14$OR),]
file14$num=1:dim(file14)[1]
p4=ggplot(file14,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+
  scale_y_continuous(breaks = file14$num,labels = file14$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p4

tmp_overlap=intersect(file11$TFname,file12$TFname)
tmp_overlap=intersect(tmp_overlap,file13$TFname)
tmp_overlap=intersect(tmp_overlap,file14$TFname)


file111=file11[file11$TFname %in% tmp_overlap,]
file111=file111[order(file111$OR),]
file111$num=1:dim(file111)[1]

p11=ggplot(file111,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[1])+
  scale_y_continuous(breaks = file111$num,labels = file111$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())

anno=data.frame(TFname=file111$TFname,num=file111$num)

file122=file12[file12$TFname %in% tmp_overlap,]
file122=file122[,-11]
file122=merge(file122,anno,by="TFname")

p12=ggplot(file122,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[2])+
  scale_y_continuous(breaks = file122$num,labels = file122$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p12


file133=file13[file13$TFname %in% tmp_overlap,]
file133=file133[,-11]
file133=merge(file133,anno,by="TFname")

p13=ggplot(file133,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[3])+
  scale_y_continuous(breaks = file133$num,labels = file133$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p13


file144=file14[file14$TFname %in% tmp_overlap,]
file144=file144[,-11]
file144=merge(file144,anno,by="TFname")

p14=ggplot(file144,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[4])+
  scale_y_continuous(breaks = file144$num,labels = file144$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p14

library(ggpubr)
ggarrange(p11,p12,p13,p14,ncol=2,nrow=2)

file2$`1/OR`=1/file2$OR
file21=file2[file2$con.file==cf[1],]
file22=file2[file2$con.file==cf[2],]
file23=file2[file2$con.file==cf[3],]
file24=file2[file2$con.file==cf[4],]

tmp_overlap2=intersect(file21$TFname,file22$TFname)
tmp_overlap2=intersect(tmp_overlap2,file23$TFname)
tmp_overlap2=intersect(tmp_overlap2,file24$TFname)

file211=file21[file21$TFname %in% tmp_overlap2,]
#file211=file211[order(file211$OR),]
file211=file211[order(file211$p.value,decreasing = T),]
file211$num=1:dim(file211)[1]

p21=ggplot(file211,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[1])+
  scale_y_continuous(breaks = file211$num,labels = file211$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())

anno2=data.frame(TFname=file211$TFname,num=file211$num)

file222=file22[file22$TFname %in% tmp_overlap2,]
file222=merge(file222,anno2,by="TFname")

p22=ggplot(file222,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[2])+
  scale_y_continuous(breaks = file222$num,labels = file222$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p22


file233=file23[file23$TFname %in% tmp_overlap2,]
file233=merge(file233,anno2,by="TFname")

p23=ggplot(file233,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[3])+
  scale_y_continuous(breaks = file233$num,labels = file233$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p23


file244=file24[file24$TFname %in% tmp_overlap2,]
file244=merge(file244,anno2,by="TFname")

p24=ggplot(file244,aes(`1/OR`,num))+geom_point(aes(size=counts,color=p.value))+
  scale_color_gradient(low = "red", high = "green")+labs(x="1/OR",y="TF.name",title = cf[4])+
  scale_y_continuous(breaks = file244$num,labels = file244$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
p24

library(ggpubr)
ggarrange(p21,p22,p23,p24,ncol=2,nrow=2)

########################文章用图
library(ggrepel)
file=read.csv("enrichment_by_TF_predit_by_MBR.csv",head=T)
file2=file[file$FDR<0.1&file$OR>1,]
file2=file2[file2$con.file=="nobias_AShM",]
file2$`1/OR`=1/file2$OR
#file2$num=as.numeric(c(1:7,9:15,8,16:20))
file2$num=as.numeric(c(4:8,1:2,9:15,3,16:20))
file2=file2[order(file2$num),]
file2$num=20:1
data1=data.frame(TFname=file2$TFname,num=file2$num)
data1$`1/OR`=file2$`1/OR`
 data2=data1
 data2$`1/OR`=0
 data_line=rbind(data1,data2)
 
ggplot(file2,aes(`1/OR`,num))+geom_point(aes(size=motif_AShM,color=FDR))+geom_text_repel(aes(label = motif_AShM))+
  scale_color_gradient(low = "red", high = "green")+geom_line(data = data_line,aes(x=`1/OR`,y=num,group=TFname))+labs(x="1/OR",y="TF.name")+
  scale_y_continuous(breaks = file2$num,labels = file2$TFname)+
  theme_bw()+theme(panel.grid.minor = element_blank())
  
  layout=matrix(c(1,1,2,2,1,1,3,3,1,1,4,4),3,4,byrow=T)
  multiplot(plotlist=list(p1,p2_22,p3,p4),layout=layout)
  
  
